#ifndef _core_math_basic_matrix_hpp_
#define _core_math_basic_matrix_hpp_

/** \file basic_matrix.hpp Contains the template class basic_matrix, that represents a mathematical basic_matrix.
 *  \author Adrian Schweizer
 *  \created  $Mo 23 Aug 10:32:12 pm CNST 2007 schwadri@SchwadriComp.local$
 *  \modified $So 28 Okt 04:32:20 pm CET 2007 schwadri@SchwadriLaptop$
 *  \todo
 *  -   Support for compile time array format selection (c or fortran layout)
 */

#include <algorithm>
#include <cmath>

#include <boost/preprocessor/control/expr_if.hpp>
#include <boost/preprocessor/empty.hpp>
#include <boost/preprocessor/identity.hpp>
#include <boost/preprocessor/cat.hpp>
#include <boost/preprocessor/repetition.hpp>
#include <boost/preprocessor/arithmetic/add.hpp>
#include <boost/preprocessor/punctuation/comma_if.hpp>
#include <boost/preprocessor/iteration/local.hpp>
#include <boost/preprocessor/facilities/expand.hpp>
#include <boost/preprocessor/comparison/greater_equal.hpp>
#include <boost/preprocessor/comparison/not_equal.hpp>
#include <boost/preprocessor/logical/and.hpp>
#include <boost/preprocessor/tuple/elem.hpp>
#include <boost/serialization/access.hpp>
//#include <boost/serialization/nvp.hpp>

#include "init.hpp"

#ifndef CORE_MATH_LIMIT_MATRIX_DIMENSION
#define CORE_MATH_LIMIT_MATRIX_DIMENSION   4
#endif

namespace core {

    namespace math {

        //forward decl
        template <class T, int M, int N=M>
            class matrix;

        template <class T, int M>
            class vector;


        template<class T, int M, int N>
            struct basic_matrixRowProxy
            {
                public:
                    basic_matrixRowProxy(T* data)
                        :   m_data(data)
                    { }
                    T& operator [](int i)
                    {
                        return *(m_data+i*M);
                    }
                private:
                    T*    m_data;
            };

        template<class T, int M, int N>
            struct Constbasic_matrixRowProxy
            {
                public:
                    Constbasic_matrixRowProxy(const T* data)
                        :   m_data(data)
                    { }
                    T operator [](int i) const
                    {
                        return *(m_data+i*M);
                    }
                private:
                    const T*    m_data;
            };

        namespace format {

            template<int M, int N>
                struct c
                {
                    template<class Index>
                        int index_to_pos(const Index& i)
                        {
                            return N*i[0] + i[1];
                        }
                };

            template<int M, int N>
                struct fortran
                {
                    template<class Index>
                        int index_to_pos(const Index& i)
                        {
                            return i[0] + M*i[1];
                        }
                };

            template<int M, int N>
                struct opengl : public fortran<M,N> { };
        }

        template<class T, int M, int N, class Derived>
            class basic_matrix
            {
            public:
                typedef T                               value_type;
                typedef basic_matrix<T,M,N, Derived>    my_type;
                typedef value_type*                     pointer;
                typedef const value_type*               const_pointer;
                typedef pointer                         iterator;
                typedef const_pointer                   const_iterator;
                typedef value_type&                     reference;
                typedef const value_type&               const_reference;

                basic_matrix()
                { }

                ~basic_matrix()
                { }

                explicit basic_matrix(const value_type& value)
                {
                    using std::uninitialized_fill;
                    uninitialized_fill(m_data,m_data+M*N,value);
                }

                /*basic_matrix(const Derived& o)
                {
                    using std::copy;
                    copy(o.begin(),o.end(),begin());
                }*/

                template<class iter>
                    basic_matrix(iter front)
                    {
                        using std::copy;
                        copy(front,front+M*N,begin());
                    }


                /*Derived&    operator=(const my_type& m)
                {
                    using std::copy;
                    copy(m.begin(),m.end(),begin());
                    return static_cast<Derived&>(*this);
                }*/

                const_iterator  begin() const
                {
                    return m_data;
                }

                inline const_iterator  end() const
                {
                    return &m_data[M*N];
                }

                iterator  begin()
                {
                    return m_data;
                }
                iterator  end()
                {
                    return &m_data[M*N];
                }

                basic_matrixRowProxy<value_type,M,N> operator[](int i)
                {
                    return basic_matrixRowProxy<value_type,M,N>(m_data+i);
                }

                Constbasic_matrixRowProxy<value_type,M,N>  operator[](int i) const
                {
                    return Constbasic_matrixRowProxy<value_type,M,N>(m_data+i);
                }


                template<class Index>
                    value_type& operator ()(const Index& i)
                    {
                        return *(m_data+i[0]+i[1]*M);
                    }

                template<class Index>
                    const value_type& operator ()(const Index& i) const
                    {
                        return *(m_data+i[0]+i[1]*M);
                    }

                template <class S>
                    Derived&  operator *= (const S& scalar)
                    {
                        for(iterator aij = begin();aij!=end();++aij)
                            *aij    *=  scalar;
                        return static_cast<Derived&>(*this);
                    }

                template <class S>
                    Derived&  operator += (const S& scalar)
                    {
                        for(iterator aij = begin();aij!=end();++aij)
                            *aij    +=  scalar;
                        return static_cast<Derived&>(*this);
                    }

                template <class S>
                    Derived&  operator -= (const S& scalar)
                    {
                        for(iterator aij = begin();aij!=end();++aij)
                            *aij    -=  scalar;
                        return static_cast<Derived&>(*this);
                    }

                template <class S>
                    Derived&  operator /= (const S& scalar)
                    {
                        for(iterator aij = begin();aij!=end();++aij)
                            *aij    /=  scalar;
                        return static_cast<Derived&>(*this);
                    }

                    Derived&  operator += (const Derived& m)
                    {
                        const_iterator bij = m.begin();
                        for(iterator aij = begin();aij<end();++aij, ++bij)
                            *aij    +=  static_cast<value_type>(*bij);
                        return static_cast<Derived&>(*this);
                    }

                Derived&  operator -= (const Derived& m)
                {
                    const_iterator bij = m.begin();
                    for(iterator aij = begin();aij<end();++aij, ++bij)
                        *aij    -=  *bij;
                    return static_cast<Derived&>(*this);
                }


                Derived& operator *= (const matrix<value_type,N,N>&  right) //2 Matrizen miteinander multiplizieren
                {
                    for(int i=0;i<M;++i)
                    {
                        value_type  temp_row[M];
                        for(int j=0;j<N;++j)
                        {
                            temp_row[j]=0.0f;
                            for(int register u=0;u<N;++u)
                                temp_row[j] += m_data[i+u*M]*right.m_data[u+j*N];
                        }
                        for(int register e=0;e<N;++e)
                            m_data[i+e*N]   =   temp_row[e];
                    }
                    return static_cast<Derived&>(*this);
                }

                vector<value_type,M> operator *(const vector<value_type,N>& v) const
                {
                    vector<value_type,M> result;
                    for(register int i1=0;i1<M;++i1)
                    {
                        result[i1]=0;
                        for(int i2=0;i2<N;++i2)
                            result[i1]+= m_data[i1+i2*M]*v[i2];
                    }
                    return result;
                }
                //static methods for creation

            private:

                //serialization for boost::archive
                friend class boost::serialization::access;
                template<class Archive>
                    void serialize(Archive & ar, const unsigned int version)
                    {
                        //ar & BOOST_SNRIALIZATION_NVP(m_data);
                        for(iterator iter=begin();iter!=end();++iter)
                            ar & *iter;
                    }

                friend struct set_entry<Derived,1,M*N>;

                template<class S>
                    void set(const S& value)
                    {
                        using std::fill;
                        fill(m_data,&m_data[M*N],value);
                    }
            protected:
                value_type       m_data[M*N];
            };

        template<class T, int M, int N>
            class matrix;

        #define TEMPLATE_PARAM_INNER(z,n2,n1)       class BOOST_PP_CAT(T,BOOST_PP_CAT(n1,n2))
        #define TEMPLATE_PARAM_ADDER(z,n1,n2)       BOOST_PP_ENUM( n2, TEMPLATE_PARAM_INNER, n1 )

        #define PARAM_ARG_INNER(z,n1,n2)            BOOST_PP_CAT(T,BOOST_PP_CAT(n2,n1)) BOOST_PP_CAT(val,BOOST_PP_CAT(n2,n1))
        #define GENERATE_PARAM_ARG_ADDER(z,n1,n2)   BOOST_PP_ENUM( n2, PARAM_ARG_INNER, n1 )

        #define PARAM_ARG2_INNER(z,n2,n1)           this->m_data[ n2 +BOOST_PP_TUPLE_ELEM(2,1,n1)*BOOST_PP_TUPLE_ELEM(2,0,n1)] = static_cast<T>(BOOST_PP_CAT(val,BOOST_PP_CAT(n2,BOOST_PP_TUPLE_ELEM(2,1,n1))));
        #define GENERATE_PARAM_ARG2_ADDER(z,n1,n2)  BOOST_PP_REPEAT( n2, PARAM_ARG2_INNER, (n2,n1) )

        #define EAT_COOKIE(n)
        #define CORE_GENERATE_OBJECT_TEMPLATE_SPEC(n)        BOOST_PP_IF(BOOST_PP_AND(BOOST_PP_GREATER_EQUAL(n,2),BOOST_PP_NOT_EQUAL(n,4)),CORE_GENERATE_OBJECT_TEMPLATE_SPEC2,EAT_COOKIE) (n)

        #define CORE_GENERATE_OBJECT_BASE_SPEC(n)                                                       \
                typedef basic_matrix<T,n,n,matrix<T,n> > base_type;                                     \
				typedef matrix<T,n,n>					 type;											\
            public:                                                                                     \
                matrix() {}                                                                             \
                                                                                                        \
                matrix(const T& value)                                                                  \
                :   base_type(value)                                                                    \
                { }                                                                                     \
                                                                                                        \
                template<BOOST_PP_ENUM(n,TEMPLATE_PARAM_ADDER,n)>                                       \
                    matrix(BOOST_PP_ENUM(n,GENERATE_PARAM_ARG_ADDER,n))                                 \
                    {                                                                                   \
                        BOOST_PP_REPEAT(n,GENERATE_PARAM_ARG2_ADDER,n)                                  \
                    }                                                                                   \
																										\
				template<class S>																		\
                    set_entry<type,1,n*n> operator=(const S& value)										\
                    {																					\
                        this->m_data[0] = value;														\
                        return set_entry<type,1,n*n>(*this);											\
                    }																					

        #define CORE_GENERATE_OBJECT_TEMPLATE_SPEC2(n)                                                  \
            template<class T>                                                                           \
                class matrix<T,n,n> : public basic_matrix<T,n,n, matrix<T,n,n> >                        \
                {                                                                                       \
                    CORE_GENERATE_OBJECT_BASE_SPEC(n)                                                   \
                };

        #define BOOST_PP_LOCAL_MACRO   CORE_GENERATE_OBJECT_TEMPLATE_SPEC
        #define BOOST_PP_LOCAL_LIMITS (2, CORE_MATH_LIMIT_MATRIX_DIMENSION)
        //DO NOT REMOVE WHITESPACE   ^
        #include BOOST_PP_LOCAL_ITERATE()


        //partial class specialization for 4x4 Matrizes

        template<class T>
            class matrix<T,4,4> : public basic_matrix<T,4,4, matrix<T,4,4> >
            {
                CORE_GENERATE_OBJECT_BASE_SPEC(4)

                //helper for homogenous transformations. transform a 2 M vector with a 3x3 basic_matrix, or a 3 M vector with a 4x4 basic_matrix
                /*vector<value_type,M-1> operator *(vector<value_type,M-1>& v)
                  {
                      vector<value_type,M-1> result;
                      for(register int i1=0;i1<M;++i1)
                      {
                      result[i1]=0;
                      for(register int i2=0;i2<(M-1);++i2)
                      {
                      result[i1]+= m_data[i1+i2*M]*v[i2];
                      }
                      result[i1]+=m_data[i1+(M-1)*M];
                      }
                      return result;
                  }
                */

                //static methods for creation

                static matrix<T,4,4> identity()
                {
                    return matrix<T,4,4>( 1, 0, 0, 0,
                                        0, 1, 0, 0,
                                        0, 0, 1, 0,
                                        0, 0, 0, 1);
                }

                template<class S>
                    static matrix<T,4,4> rotatex(const S& angle)
                    {
                        using std::cos;
                        using std::sin;
                        T c = static_cast<T>(cos(angle));
                        T s = static_cast<T>(sin(angle));

                        return matrix<T,4,4>( 1, 0, 0, 0,
                                            0, c,-s, 0,
                                            0, s, c, 0,
                                            0, 0, 0, 1);
                    }

                template<class S>
                    static matrix<T,4,4> rotatey(const S& angle)
                    {
                        using std::cos;
                        using std::sin;
                        T c = static_cast<T>(cos(angle));
                        T s = static_cast<T>(sin(angle));

                        return matrix<T,4,4>( c, 0, s, 0,
                                            0, 1, 0, 0,
                                           -s, 0, c, 0,
                                            0, 0, 0, 1);
                    }

                template<class S>
                    static matrix<T,4,4> rotatez(const S& angle)
                    {
                        using std::cos;
                        using std::sin;
                        T c = static_cast<T>(cos(angle));
                        T s = static_cast<T>(sin(angle));

                        return matrix<T,4,4>( c,-s, 0, 0,
                                            s, c, 0, 0,
                                            0, 0, 1, 0,
                                            0, 0, 0, 1);
                    }

                static matrix<T,4,4> mirrorx()
                {
                    return matrix<T,4,4>(-1, 0, 0, 0,
                                        0, 1, 0, 0,
                                        0, 0, 1, 0,
                                        0, 0, 0, 1);
                }

                static matrix<T,4,4> mirrory()
                {
                    return matrix<T,4,4>( 1, 0, 0, 0,
                                        0,-1, 0, 0,
                                        0, 0, 1, 0,
                                        0, 0, 0, 1);
                }

                static matrix<T,4,4> mirrorz()
                {
                    return matrix<T,4,4>( 1, 0, 0, 0,
                                        0, 1, 0, 0,
                                        0, 0,-1, 0,
                                        0, 0, 0, 1);
                }


                static matrix<T,4,4> scale(const vector<T,3>& scale)
                {
                    return matrix<T,4,4>(scale[0], 0, 0, 0,
                                        0, scale[1], 0, 0,
                                        0, 0, scale[2], 0,
                                        0,      0,   0, 1);
                }

                template<class S>
                    static matrix<T,4,4> scale(const S& scale)
                    {
                        return matrix<T,4,4>(scale, 0, 0, 0,
                                            0, scale, 0, 0,
                                            0, 0, scale, 0,
                                            0, 0,   0,   1);
                    }

                static matrix<T,4,4> translate(const vector<T,3>& v)
                {
                    return matrix<T,4,4>( 1, 0, 0, v[0],
                                        0, 1, 0, v[1],
                                        0, 0, 1, v[2],
                                        0, 0, 0, 1);
                }

                template<class S>
                    static matrix<T,4,4> frustum(const vector<T,2>& bl,const vector<T,2>& ur,const S& near, const S& far)
                    {
                        return matrix<T,4,4>(2*near/(ur[0]-bl[0]), 0, (ur[0]+bl[0])/(ur[0]-bl[0]), 0,
                                            0, 2*near/(ur[1]-bl[1]), (ur[1]+bl[1])/(ur[1]-bl[1]),0,
                                            0, 0, -(far+near)/(far-near), -2*far*near/(far-near),
                                            0, 0, 1, 0);
                    }

                template<class S>
                    static matrix<T,4,4> frustumsym(const vector<T,2>& size,const S& near, const S& far)
                    {
                        return matrix<T,4,4>(2*near/(size[0]), 0, 0, 0,
                                            0, 2*near/(size[1]), 0, 0,
                                            0, 0, -(far+near)/(far-near), -2*far*near/(far-near),
                                            0, 0, 1, 0);
                    }

                template<class S>
                    static matrix<T,4,4> ortho(const vector<T,2>& bl,const vector<T,2>& ur,const S&  near, const S& far)
                    {
                        return matrix<T,4,4>(2/(ur[0]-bl[0]), 0, 0, 0,
                                            0, 2/(ur[1]-bl[1]), 0, 0,
                                            0, 0, -2/(far-near), -2*(far+near)/(far-near),
                                            0, 0, 0, 1);
                    }

                template<class S>
                    static matrix<T,4,4> perspective(const S& fov, const S& aspect, const S& near, const S&  far)
                    {
                        return matrix<T,4,4>( 1.0/(tan(fov/2.0)*aspect), 0, 0, 0,
                                                0, 1.0/(tan(fov/2.0)), 0, 0,
                                                0, 0, -(far+near)/(far-near), -2*far*near/(far-near),
                                                0, 0, 1, 0);
                    }
            };

        #undef CORE_GENERATE_OBJECT_BASE_SPEC
        #undef CORE_GENERATE_OBJECT_TEMPLATE_SPEC
        #undef CORE_GENERATE_OBJECT_TEMPLATE_SPEC2

        #undef TEMPLATE_PARAM_ADDER
        #undef TEMPLATE_PARAM_INNER
        #undef PARAM_ARG_INNER
        #undef PARAM_ARG2_INNER
        #undef GENERATE_PARAM_ARG_ADDER
        #undef GENERATE_PARAM_ARG2_ADDER
        #undef EAT_COOKIE
        //usefull shortcuts

        typedef matrix<float,4,4> matrix4f;
        typedef matrix<float,3,3> matrix3f;
        typedef matrix<float,2,2> matrix2f;

        typedef matrix<double,4,4> matrix4d;
        typedef matrix<double,3,3> matrix3d;
        typedef matrix<double,2,2> matrix2d;

        //serialization methods
        template<class CharT, class CharTraits, class T, int M, int N>
            std::basic_ostream<CharT,CharTraits>& operator << (std::basic_ostream<CharT,CharTraits>& out, const matrix<T,M,N>& m)
            {
                //FIXME: this only works with char

                out << '[' << *m.begin();

                for(typename matrix<T,M,N>::const_iterator miter=m.begin()+1;miter!=m.end();++miter)
                    out << ',' << *miter;

                out << ']';
                //using std::copy;
                //copy(v.begin(),v.end(),std::ostream_iterator<T,CharT,CharTraits>(out,","));
                return out;
            }

    } // namespace math

} // namespace core

#include "matrix_ops.hpp"
#endif // _basic_matrix_hpp_
